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Abstract 

Background: There has been a considerable increase in studies investigating rates of diversification and character 
evolution, with one of the promising techniques being the BiSSE method (binary state speciation and extinction). 
This study uses simulations under a variety of different sample sizes (number of tips) and asymmetries of rate 
(speciation, extinction, character change) to determine BiSSE's ability to test hypotheses, and investigate whether 
the method is susceptible to confounding effects. 

Results: We found that the power of the BiSSE method is severely affected by both sample size and high tip ratio 
bias (one character state dominates among observed tips). Sample size and high tip ratio bias also reduced 
accuracy and precision of parameter estimation, and resulted in the inability to infer which rate asymmetry caused 
the excess of a character state. In low tip ratio bias scenarios with appropriate tip sample size, BiSSE accurately 
estimated the rate asymmetry causing character state excess, avoiding the issue of confounding effects. 

Conclusions: Based on our findings, we recommend that future studies utilizing BiSSE that have fewer than 
300 terminals and/or have datasets where high tip ratio bias is observed (i.e., fewer than 10% of species are 
of one character state) should be extremely cautious with the interpretation of hypothesis testing results. 

Keywords: Key innovations, Character evolution, Systematics 



Background 

Maddison et al. [1] described a method using a binary- 
state speciation and extinction model (BiSSE) that esti- 
mates rates of change in a binary character and rates of 
speciation and extinction contingent on the character 
state, given a known distribution of observed states on 
the tips of a tree of contemporaneous species. The BiSSE 
method assumes that the tree includes all extant species 
and that all species have known data for the state of the 
single binary character [1]. Fitzjohn et al. [2] provided 
additional methodology for including known species di- 
versity into an incomplete phylogeny if a researcher 
could confidently place taxa into unresolved terminal 
clades. BiSSE provides estimates for the rates of speci- 
ation in each character state (A 0 > ^i)> extinction in each 
state (u 0 , (ii), and character transition rates between 
states (qoi> <7io)- Such estimates are important to studies 
of whether a particular feature is controlling the diversi- 
fication rates of clades and whether the effect is on 
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speciation, extinction, or both [3-5]. Recent studies have 
also modified the basic approach of the BiSSE method to 
estimate further parameters associated with quantitative 
characters [6] and geography [7]. 

Maddison et al. [1] discussed the development of the 
BiSSE method and demonstrated its ability to estimate 
rates from simulated data sets. Recently there has been a 
burst in the number of studies that have utilized the 
BiSSE method to explore rates of diversification in vari- 
ous taxonomic groups for the purpose of testing hypoth- 
eses involving key innovations and the evolution of 
characters [8-19]. However, the majority of these studies 
explore diversification and character evolution hypoth- 
eses with fewer than 200 taxa, despite the initial warning 
by Maddison et al. [1] that the power of analysis may be 
affected by low sample size. 

Wise use of any statistical method should be guided 
by an understanding of its power and ability to distin- 
guish hypotheses of interest, but current empirical 
studies lack sufficient guidance, because there has been 
little work on the BiSSE method's behavior under dif- 
ferent sample sizes (numbers of species) and parameter 
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values. The primary goal of this study is to explore the 
power and accuracy of parameter estimation of the 
BiSSE method. Using simulations we explore the num- 
ber of species needed to obtain good power, the 
advantages of estimating fewer parameters, and the ef- 
fect of the extreme asymmetries in rates. Additionally, 
because there are many ways that an observed excess 
of a character state can be explained through macro- 
evolutionary processes (e.g., increased speciation rates 
in taxa with State 1, higher extinction rates in taxa 
with State 0), there is also concern regarding con- 
founding effects [20], and whether BiSSE can identify 
which rate asymmetries are causing the observed char- 
acter excess in empirical data. 

Results 

Power of BiSSE method 
Asymmetries in speciation rate 

Extremely low power (<5%) was observed when tree size 
included 50 taxa, regardless of the degree of rate asym- 
metry (Figure la, Additional file 1: Table SI). Power 
marginally improved for tree sizes of 100 taxa, but 
decreased considerably as the asymmetry increased to 
10 and 20x rate difference (Figure la). A tree size of 300 
tips indicated a higher overall power for each difference 
in rate than those observed with 50 or 100 taxa. Power 
increased as the degree of difference in rate asymmetry 
increased, until reaching a rate difference of four times 
the speciation rate where the power begins to decrease 
as the degree of rate asymmetry grows. This same pat- 
tern was observed for a tree size of 500 tips (Figure la). 
In general, power increased as tree size increased 
(Figure la), and a pattern of power decrease following 
an increase in tip-ratio bias resulting from rate asym- 
metry was observed in the simulation including more 
tips (Figure la). 

Asymmetries in rate of character state change 

For each asymmetrical model of character change (e.g., 
2x, 5x), power increased with an increase in tree size 
(Figure lb). Power remained low and did not increase 
with a difference in rates when the tree size was 50 taxa. 
There was a slight increase in power as the degree of dif- 
ference in rates increased with tree sizes of 100 taxa 
(Figure lb). Power was higher for simulations with 300 
taxa for each respective difference in rates compared to 
the same simulations with 100 and 50 taxa. Power 
increased as the rate difference increased to 5x then lev- 
eled off to lOx, followed by a strong decrease in power 
as the difference in rates increased to 20x and 40x 
(Figure lb). This same pattern was observed in simula- 
tions of 500 taxa and in general there was a decrease in 
power observed in all tip sizes beyond a lOx difference 
in rates of character change (Figure lb). 
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Figure 1 Power of BiSSE under simulations with asymmetrical 
rates of (a) speciation, (b) character change, (c) extinction. See 

Table 1 and Additional file 1: Table S1-S3 for list of rate values. 



Asymmetries in extinction rate 

As with rates of speciation and character state change, 
power increased as tree size increased regardless of the 
amount of difference in extinction rate (Figure lc). With 
tree sizes of 50 taxa, power was low regardless of the de- 
gree of rate asymmetry, and power was similarly low 
with 100 taxa. With tree sizes of 300 and 500 taxa, 
power increased until a rate difference of 3x, with power 
decreasing as the degree of rate asymmetry continue to 
increase. Overall, power in hypothesis testing with ex- 
tinction rates was lower than those for speciation or 
character state change (Figure lc). 
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Power of six parameter model vs. four parameter model 

The change in power among low and high tip bias sce- 
narios when using a model that estimates fewer para- 
meters is observed in Figure 2 (Additional file 1: Table 
S4). The reduced four parameter model differs from the 
full six parameter model (A. 0) \ 1( u 0 , Hi, q 01 , q 10 ) with two 
of the rates constrained to be equal in both character 
states, for example (\q, \i, |io=Ui, <?oi=<7io)- In nearly all 
cases, especially with tree sizes greater than 300 taxa, 
there was an increase in power when the four parameter 
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Figure 2 Power of BiSSE under four versus six parameter 
models under low (3,:1 0 ) and high (7,:1 0 ) character state bias 
as the result of rate asymmetries in (a) speciation, (b) character 
change, and (c) extinction. See Additional file 1: Table S4 for 
rate values. 



model was used compared to the six parameter model 
(Figure 2). The greatest increase in power occurred in the 
scenarios with a greater degree of bias (7:1 State 1 favored), 
and the effect was observed regardless of which parameter 
possessed the rate asymmetry (Figure 2, Additional file 1: 
Table S4). As was the case with the six parameter estima- 
tions, power is low when tree size is below 300 taxa with a 
four parameter model. 

Parameter estimation 

Estimating parameters under asymmetrical speciation 

As previously identified by Maddison et al. (2007), the 
BiSSE method estimates speciation rates well, with 
strong delimitation of known asymmetrical rates given 
an appropriate tree size (Figure 3a), although precision 
decreases as the speciation rate asymmetry, and tip bias, 
increases. For estimates of known symmetrical rates of 
character change under asymmetrical speciation rates, 
accuracy and precision of estimating the known rates be- 
come significantly worse under the higher tip bias sce- 
nario (Figure 3b). A similar pattern is recovered for 
estimates of known symmetrical extinction values where 
rates were estimated with more accuracy and precision 
under low tip bias (asymmetrical speciation rate of 1.25), 
with accuracy and precision strongly decreasing under 
high tip bias (20 x asymmetry in speciation rates) as 
observed in Figure 3c. Accuracy and precision of param- 
eter estimation greatly decreases for all rates with a 
small sample size of tips, regardless of low or high tip 
bias (Additional file 2: Figure SI). 

Estimating parameters under asymmetrical character change 

Estimates of asymmetries in character change are not as 
accurate or precise as estimating rates of speciation 
(Figure 4b). In general, precision decreases as the rate 
difference increase results in a high tip bias scenario, 
with parameter estimation being poor for a known esti- 
mate of a 40x rate difference in character change. Sym- 
metrical speciation rates (\ 0 = 0.1, \±= 0.1) are well 
estimated when the rate of character change is 2x (low 
tip bias) as seen in Figure 4a. However, with a 40 x (high 
tip bias) difference in the rate of character change the 
precision of parameter estimation appears to decrease, 
and the number of estimates for highly asymmetrical 
speciation rates increases (Figure 4a). Parameter esti- 
mation of known extinction rates (^n = 0.03, fii = 
0.03) is more accurate under the 2x (low tip bias) sce- 
nario rather than the 40x (high tip bias) scenario 
where estimates of known extinction values are very 
poor (Figure 4c). 

Estimating parameters under asymmetrical extinction 

Estimates of known extinction values are poor and 
seem to lack precision, with precision decreasing as 
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Figure 3 Parameter estimations of (a) speciation, (b) character change, and (c) extinction under different degrees of asymmetry in 
speciation rates with corresponding tip ratios. Point of intersection between red lines represents known values. 



the difference in extinction rates increase (Figure 5c). 
Speciation values are well estimated to the known 
values (Xo = 0.1, Xi= 0.1) when the difference between 
extinction rates is 2x (low tip bias), but accuracy and 
precision seems to decrease dramatically as the extinc- 
tion rate asymmetry increases (high tip bias), leading 
to an abundance of highly asymmetrical estimates of 
speciation rates (Figure 5a). Parameter estimation of 
character change has the same pattern, in which the 
known rate {q 01 = 0.01, q 10 = 0.01) is well estimated 
under 2x (low tip bias) extinction rate differences 
(Figure 5b), but is poorly estimated under an increased 
difference in extinction rate asymmetry that leads to 
high tip bias (Figure 5b). This poor estimation leads to 
a dramatic increase in asymmetrical rates of character 
change that favor transitions from State 0 to 1. 



Discussion 

Impact of tree size on power 

The statistical power of the BiSSE method depends on 
the number of taxa and the degree of asymmetry in rates 
of speciation, extinction, and character state change. In 
terms of tree size, BiSSE achieves extremely low power 
when testing hypotheses of rate asymmetry if fewer than 
100 taxa are used in the analysis, even when rates are 
known to be highly asymmetrical (Figure 1). As a result, 
the potential for a Type II error (failing to reject the null 
hypothesis when the alternate hypothesis is true) is ex- 
tremely high. 

The highest power attributed to any rate asymmetry 
associated with 300 taxa is only 50% (Figure 1) under a 
six parameter model. Researchers that attempt to utilize 
the BiSSE method with fewer than 300 taxa should take 
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Figure 4 Parameter estimations of (a) speciation, (b) character change, and (c) extinction under different degrees of asymmetry in 
character change rates with corresponding tip ratios. Point of intersection between red lines represents known values. 



caution. Below 200 taxa there is little statistical power 
associated with identifying rate asymmetries, regardless 
of whether strong asymmetries exist. Maddison et al. [1] 
hypothesized that large amounts of data would be 
needed to distinguish significant asymmetries because 
there are many ways to arrive at a given phylogeny. 

Power of simplified models 

If external information justifies a simplification of the 
model, then greater power can be achieved. We found an 
overall increase in power when utilizing a four parameter 
model over a six parameter model regardless of which 
rate possesses the asymmetry (Figure 2, Additional file 1: 
Table S4), although tree sizes greater than 300 tips are 
still desirable. Whether researchers can take advantage of 
this greater power using a model with fewer free para- 
meters depends of course on whether it is reasonable to 



assume in advance that any asymmetries are restricted to 
one rate. 

Confounding processes 

Strong asymmetries in rates of speciation, character 
state change, and extinction yielded, as expected, a 
strong excess of tips with a single character state 
(Table 1). The question to a biologist observing such a 
pattern is; which rate's asymmetry might have led to 
such excess? Maddison [20] hypothesized that teasing 
apart parameters that are the cause of taxonomic state 
frequency asymmetry is difficult, and that simultaneously 
estimating these parameters may help address this issue. 
We found that under best-case scenarios where high tip 
bias is absent and sample size is high (preferably greater 
than 300 taxa), parameters are estimated accurately and 
precisely (Figures 3, 4, 5, Additional file 2: Figure SI, 
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Figure 5 Parameter estimations of (a) speciation, (b) character change, and (c) extinction under different degrees of asymmetry in 
extinction rates with corresponding tip ratios. Point of intersection between red lines represents known values. 



Table 2). This indicates that the BiSSE method is cap- 
able of identifying the rate asymmetry that is causing 
taxonomic excess when datasets are used with a suffi- 
cient amount of tips and where high tip bias is absent, 
as seen in Figures 3, 4, 5 (Table 2). In these scenarios, 
BiSSE is able to properly identify the known process 
with the rate asymmetry (Table 2), whether it was spe- 
ciation (Figure 3a), character change (Figure 4b), or ex- 
tinction (Figure 5c). This is further corroborated by the 
power, which remains high in cases where sample size 
is large and tip ratio bias is minimized (Figure 1). How- 
ever, teasing apart which rate asymmetries are causing 
taxonomic excess of one character state is a problem 
under two conditions; if tree size is small (Additional 
file 2: Figure SI), regardless of the degree of tip bias, 
and if there is a high degree of tip ratio bias towards 



one state, even while simultaneously estimating para- 
meters using the BiSSE method with trees of substantial 
size (500 tips, Figures 3, 4, 5, Additional file 2: Figure SI, 
Table 2). Below we discuss why high tip ratio bias leads 
to a decrease in power and parameter estimation, and 
how much tip ratio bias is necessary to have an adverse 
affect on teasing apart confounding processes. 

When taxa are saturated for a particular state (high tip 
ratio bias), the BiSSE method estimates high rate asym- 
metries to explain this pattern, even for those rates 
known to be low and symmetrical (Figure 3, 4, 5, 
Table 2). For example, when extinction is lOx higher for 
State 0 than State 1 (Tip bias of IO^Iq), the extinction 
asymmetry is detected, but the rate of character change 
is also estimated to be highly asymmetrical (Table 2). 
The method infers rapid change from State 0 to 1 
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Table 1 Parameter values for rate asymmetry simulations 
used to assess power and parameter estimation of BiSSE 
method 



Rate asymmetry 


Tip ratio 


Figure 


% State 0 


Speciation (q 0 i and q lu = 0.01, 
/j 0 and /j, = 0.03) 










1.25x(A 0 = 0.1, A,= 0.125) 


3,: 


lo 


1A 


29.23 


1.5x(A 0 = 0.1, A,= 0.15) 


5,: 


1o 


1A 


19.33 


2x (A 0 = 0.1, A,= 0.2) 


10, 


:1o 


1A 


9.90 


3x (A 0 = 0.1, A,= 0.3) 


20, 


:1o 


1A 


4.94 


4x (A 0 = 0.1, A,= 0.4) 


30, 


:1o 


1A 


3.19 


5x (A 0 = 0.1, A,= 0.5) 


40, 


:1o 


1A 


2.50 


10x (A 0 = 0.1, A,= 1.0) 


90, 


:1o 


1A 


1.13 


20x (A 0 = 0.1, A,= 2.0) 


180 


,:lo 


1A 


0.51 


Character Change (A 0 and A,= 0.1, 
u 0 and /j, = 0.03) 










2x (t? 0 , = 0.01, q m = 0.005) 


2,: 


lo 


IB 


33.96 


3x (q 0] = 0.015, q ]0 = 0.005) 


3,: 


lo 


1B 


24.01 


4x (q m = 0.02, q 10 = 0.005) 


4,: 


lo 


IB 


20.57 


5x (q 0] = 0.025, q ]0 = 0.005) 


5,: 


lo 


IB 


16.69 


10x (g 0 i = 0-05, q 10 = 0.005) 


10, 


:1o 


1B 


9.14 


20x(q 0 i = 0.1, q, 0 = 0.005) 


20, 


:1o 


IB 


4.69 


40x (q m = 0.2, q m = 0.005) 


40, 


:1o 


1B 


2.39 


Extinction (A 0 and Ai= 0.1, q 0 , 
and q 10 = 0.01) 










2x (jj 0 = 0.06, = 0.03) 


3,: 


lo 


1C 


23.85 


3x (u c = 0.09, = 0.03) 


6,: 


lo 


1C 


13.21 


4x = 0.12, (J, =0.03) 


9,: 


lo 


1C 


9.29 


5x (fj 0 = 0.15, =0.03) 


12, 


:1o 


1C 


7.12 


1 0x (u 0 = 0.3, ^ = 0.03) 


27, 


:1o 


1C 


3.40 



The observed percent of terminals with State 0 are for simulations with 500 
taxa. Observed tip ratios were nearly identical for other taxa sizes (50, 100, 
and 300). 



(Figure 5b, Table 2) when in fact, the constrained rates 
were fairly low and symmetrical (qoi = 0.01, qio = 0.01). 
The inability to identify confounding processes when 
there is high tip ratio bias results in the observed de- 
crease in power when the degree of rate asymmetry 
increases in either the rate of speciation, extinction, or 
character change (Figure 1). This issue is not resolved 
with larger tree sizes (Figure 1), and the degree of tip 
bias needed to adversely affect the accuracy and preci- 
sion of parameter estimation varied with the process 
(e.g., speciation, extinction, character change). 

The amount of bias needed to reduce accuracy and 
precision of parameter estimation, leading to worst-case 
scenarios, changed depending on the process. When es- 
timating parameters associated with asymmetries of spe- 
ciation rates, power sharply decreases when fewer than 
2.5% of the taxa have one of the binary traits (Tables 1, 
2, Figures 1, 3). Character change and extinction rate 



asymmetries are more adversely affected by trait rarity, 
with the worst-case scenarios happening when a binary 
trait occurs in fewer than 8-10% of the terminal taxa. 
As with speciation, high tip ratio bias leads to a sharp 
decrease in power regardless of sample size (Figures lb, 
lc, 4, 5; Tables 1, 2, Additional file 1: Tables S1-S3). 
Overall, caution is recommended when using the BiSSE 
method if tip ratio bias is greater than a 10:1 ratio for 
binary states in terminal taxa, as this level of trait rarity 
in one state may have a negative impact on the power of 
the analysis (regardless of sample size) and the ability to 
identify confounding processes under these worst-case 
scenarios. 

The observed pattern of a decrease in power asso- 
ciated with speciation, extinction, and character change 
when rates become increasingly asymmetrical is troub- 
ling (Figure 1). In general, investigators should be cau- 
tious using the BiSSE method when one of the binary 
characters in question is exceedingly rare in their data 
sets (less than 10%). BiSSE may mistakenly estimate the 
wrong parameter (or combination of parameters) to be 
the cause of taxonomic excess in situations where high 
tip ratio bias is observed, and increased tree size does 
not alleviate this issue. If investigators report a signifi- 
cant result, there is a possibility that the inability to 
identify confounding processes in these scenarios may 
impact the interpretation of the causes for the observed 
pattern of state asymmetries. Further work is needed to 
establish a robust methodology within the BiSSE frame- 
work for teasing apart the parameters that are directly 
contributing to character state bias under these high tip 
ratio scenarios. 

Difference in estimation accuracy among processes of 
change 

We found that the BiSSE method is far more accur- 
ate and precise with estimates of rates of speciation 
and character change than with extinction rates 
(Table 2, Figures 3, 4, 5, Additional file 2: Figure SI), 
which was also noted by Maddison et al. [1]. Overall, 
a greater degree of tip ratio bias is needed to reduce 
the accuracy and precision of speciation parameter 
estimation and a loss of power when the rate asym- 
metry is related to speciation (discussed previously). 
While extinction can leave a signal in molecular 
phylogenetic trees recovered from extant taxa alone 
[21,22], estimating extinction rates from molecular 
phylogenies often results in rates approaching zero, 
which is potentially the result of cladogenetic events 
being directly inferred from molecular phylogenies 
and not extinction events [22]. Recently, Rabosky 
[23] indicated that without information from the fos- 
sil record, estimating extinction rates from molecular 
data alone may be potentially impossible if rates of 
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Table 2 Summary statistics of parameter estimates for rate asymmetry simulations with 500 taxa under low and high 
tip bias scenarios 



Rate asymmetry 






<7io 


<?oi 






Speciation (Figure 3) 














Low Tip Bias (3,:1 0 ) 














Simulated 


0.1 


0.125 


0.01 


0.01 


0.03 


0.03 


Estimated (Mean ± S e ) 


0.101 ± 0.018 


0.125 ± 0.013 


0.01 ± 0.003 


0.01 ± 0.006 


0.035 ± 0.031 


0.029 ± 0.022 


Medium Tip Bias (40,:1 0 ) 














Simulated 


0.1 


0.5 


0.01 


0.01 


0.03 


0.03 


Estimated (Mean ± S e ) 


0.107 ±0.1 19 


0.502 ± 0.033 


0.012 ± 0.007 


0.09 ± 0.295 


0.103 ± 0.33 


0.034 ± 0.046 


High Tip Bias (180,:1 0 ) 














Simulated 


0.1 


2.0 


0.01 


0.01 


0.03 


0.03 


Estimated (Mean ± S e ) 


0.377 ± 2.57 


2.048 ±0.124 


0.018 ± 0.028 


3.1 ± 7.6 


1.78 ±4.85 


0.11 ± 0.17 


Character Change (Figure 4) 














Low Tip Bias (2j:1o) 














Simulated 


0.1 


0.1 


0.005 


0.01 


0.03 


0.03 


Estimated (Mean ± S e ) 


0.103 ± 0.033 


0.099 ± 0.01 1 


0.005 ± 0.002 


0.0099 ± 0.005 


0.038 ± 0.053 


0.03 ± 0.02 


Medium Tip Bias (10-|:1 0 ) 














Simulated 


0.1 


0.1 


0.005 


0.05 


0.03 


0.03 


Estimated (Mean ± S e ) 


0.104 ± 0.032 


0.098 ± 0.009 


0.0058 ± 0.006 


0.057 ± 0.218 


0.048 ± 0.075 


0.026 ±0.014 


High Tip Bias (40,:1 0 ) 














Simulated 


0.1 


0.1 


0.005 


0.2 


0.03 


0.03 


Estimated (Mean ± S e ) 


0.142 ± 0.145 


0.098 ± 0.008 


0.008 ± 0.009 


0.193 ± 0.443 


0.281 ± 0.462 


0.023 ± 0.014 


Extinction (Figure 5) 














Low Tip Bias (2,:1 0 ) 














Simulated 


0.1 


0.1 


0.01 


0.01 


0.06 


0.03 


Estimated (Mean ± S e ) 


0.1 ± 0.02 


0.099 ± 0.01 


0.009 ± 0.002 


0.009 ± 0.006 


0.063 ± 0.031 


0.028 ± 0.016 


Medium Tip Bias (3 1 : 1 o) 














Simulated 


0.1 


0.1 


0.01 


0.01 


0.09 


0.03 


Estimated (Mean ± S e ) 


0.077 ±0.014 


0.102 ± 0.009 


0.009 ± 0.003 


0.012 ± 0.013 


0.093 ± 0.047 


0.027 ± 0.014 


High Tip Bias (10,:1 0 ) 














Simulated 


0.1 


0.1 


0.01 


0.01 


0.3 


0.03 


Estimated (Mean ± S e ) 


0.1 16 ±0.109 


0.098 ± 0.008 


0.013 ± 0.03 


0.194 ± 1.2 


0.29 ± 0.33 


0.028 ± 0.015 



diversification vary across a topology. As demon- 
strated by this study, BiSSE is capable of estimating 
rates of extinction from known values given sufficient 
data (e.g., large tree size, low tip bias), albeit with far 
less accuracy and precision than rates of speciation 
or character change. Also, if an asymmetry in extinc- 
tion rates leads to a high tip ratio bias, the accuracy 
and precision of extinction rates estimation decreases 
(Figure 5c). Finally, the accuracy and precision of es- 
timating rates is also impacted by a small sample size 
(<300 tips) regardless of tip ratio bias as indicated in 
Additional file 2: Figure SI. In addition to low power 
associated with small tree sizes, investigators should be 
cautious of a significant result if tree size is less than 
300 tips, as the inferred cause(s) for the potential 



evolutionary pattern may be misled by the issue of con- 
founding processes. 

Identifying critical values 

An additional issue uncovered in this analysis is the diffi- 
culty in finding appropriate critical values. The critical 
values from our simulations suffer from a great deal of 
variation (Additional file 1: Tables S1-S4). The nearly 
consistent difference in critical values suggests that sim- 
ply comparing the statistic to a x distribution may not 
be appropriate as suggested by Maddison et al. [1]. 
Therefore, we suggest simulating your own critical 
values, using as many replicates as possible, as we did 
here. Suitable simulators are available in both Mesquite 
and the Diversitree R-package [2]. 
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Implementation issues 

Another possibility is that the method performs better 
than these results suggest, but limitations of our likeli- 
hood optimizer limit the power of our implementation. 
Our implementation, as described in Maddison et al. 
[1] uses several techniques to overcome the limitations 
of the Brent [24] optimizer, including using multiple 
searches from randomly generated starting points, and 
the option of starting searches from values estimated 
from more constrained models. Tests (Midford unpub- 
lished) with alternative optimizers, using the Diversi- 
tree package described Fitzjohn et al. [2] indicate that 
little, if any, loss of power is actually related to the 
choice of optimizer. 

Conclusion 

The power of the BiSSE likelihood method to test hy- 
potheses of rate asymmetry is susceptible to both tree 
size and variation in parameter rates. Overall, power of 
the BiSSE method is low if the tree size is below 300 
taxa, and investigators should take special care to inves- 
tigate the power of their results if applying the BiSSE 
method to topologies with fewer than 300 tips. Power is 
increased when estimating fewer parameters, so utilizing 
a four parameter model to test hypotheses may be pre- 
ferable if appropriate. 

This study indicates that contrary to the hope 
expressed in Maddison [20], the problem of con- 
founding effects can still occur while estimating 
process parameters simultaneously if there is low 
sample size and/or high tip ratio bias. Under scenar- 
ios of large sample size (greater than 300 taxa) and 
low tip bias, the BiSSE method accurately and pre- 
cisely identifies the rate parameters contributing to 
the observed taxonomic excess. However, if diversifi- 
cation rate parameters are too asymmetrical (yielding 
a high tip ratio bias) and/or sample size is low, 
BiSSE is unable to accurately estimate rates. This in 
turn results in a dramatic decrease of power. We 
recommend that investigators must be cautious when 
interpreting their results if there is a character state 
bias among tips greater than a 10:1 ratio in favor of 
either binary state. In these worst-case situations, 
properly identifying the process responsible for taxo- 
nomic excess may be impossible regardless of the 
number of tips in the dataset. If investigators using 
data with fewer than 300 tips and/or with high tip 
bias report a significant result, there exists a possibil- 
ity that the issue of confounding effects has misled 
the identified rate cause(s) of the significant result. 
Further work is needed within the BiSSE framework 
to develop methods to better identify which parameters 
are causing the character state bias in these worst-case 
scenarios (e.g., low sample size, high tip bias). Further 



exploration of the impact of multiple rate asymmet- 
ries is also needed. However, it is clear that if multiple 
rate asymmetries are occurring that promote high tip 
ratio bias there will be difficulties with power and 
parameter estimation. 

Methods 

Hypothesis testing and the power of the BiSSE method 

The BiSSE likelihood calculation and parameter esti- 
mations were done in the Diverse package [25] of 
Mesquite 2.7 [26]. Maddison et al. [1] suggested that 
the probability of rejecting a false null hypothesis 
(power) may vary with the number of species in an 
analysis, and with the degree of rate difference among 
parameters. The initial exploration of power for the 
BiSSE method in Maddison et al. [1] focused on three 
rate asymmetric scenarios (one for speciation, charac- 
ter change, and extinction) with a tree size of 500 
tips. To explore the full range of power for the BiSSE 
method when using the six parameter model, 500 
trees were simulated under a variety of tip sizes and 
parameter combinations where an asymmetry in one 
rate parameter was introduced (e.g., \ 0 < 

When a tree and character are simulated with 
asymmetrical rates in one process — speciation, ex- 
tinction, or character state change — our question 
was whether the BiSSE method could detect this 
asymmetry and estimate the rates correctly. A biolo- 
gist facing such a question with real data may be 
interested in just one of the processes (e.g., is there 
an asymmetry in speciation rates?), and would there- 
fore face a choice: are the other processes assumed 
to be symmetrical (i.e. extinction rates u 0 = u.i) or 
not? Assuming the other processes to be symmetrical 
reduces the complexity of the models and permits 
the method to focus entirely on the process of inter- 
est. We studied the benefits of such simplifying 
assumptions as described in the next section. How- 
ever, because biologists typically would not have in- 
formation confirming the other processes to be 
symmetrical, in most of our analyses we permitted 
all three processes to be asymmetrical, thus requiring 
us to compare the null five-parameter model against 
a full six parameter model. 

Thus, for any given asymmetry simulation, the 
BiSSE likelihood score for the full six parameter 
model was calculated and compared to the likelihood 
score of a corresponding five parameter model where 
the rate parameter with an asymmetry being tested is 
constrained to be equal for both states (e.g., \ 0 = 
In addition, a null hypothesis set of simulations with 
500 trees was generated where all rate values are sym- 
metrical, and a distribution of the BiSSE likelihood 
score difference for the six and five parameter models 
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were calculated for each null hypothesis. Power was 
determined as the percentage of likelihood difference 
scores between the six and five parameter models of 
the asymmetrical simulations that were above the 5% 
cutoff value established by the corresponding likeli- 
hood difference score distributions of the rate sym- 
metrical null hypothesis simulations. We simulated 
sets of 500 trees with varying bias for each of the 
three parameter pairs as well as sets without bias and 
each rate parameter combination (Table 1) was tested 
under tree sizes of 50, 100, 300, and 500 taxa, re- 
spectively, in which the probability of the root state 
was stationary [27]. 

Many asymmetries in rates would yield biases in the 
frequency of the two character states observed in the 
contemporaneous species at the tips of the tree, with 
stronger asymmetries yielding strong biases. Insofar as 
BiSSE's behavior might differ depending on the strength 
of the bias in state frequencies, we explored examples 
with rate parameter asymmetries that would result in a 
range of observed high and low biases in character state 
distributions across the tips (Tables 1 and 2). Ratios of 
tip bias (taxa with state 0 relative to state 1) were calcu- 
lated from their corresponding biases in the model rates 
using the stationary frequency formula in Appendix 2 of 
Maddison et al. [1]. 

gx'(l-x')-x'q01 + (l-x")ql0 = 0 (1) 

Where g = \ 0 - u 0 - Xi+Hi, and x" is the stationary 
frequency of state 0. 

As expected, tip bias increased as asymmetries in the 
simulation parameters increased steadily across all tree 
sizes, and the observed asymmetry in taxa with each 
state matched expectations given the starting rate 
asymmetries (Table 1, Additional file 1: Tables SI, S2, 
and S3). 

Power in reduced (4-parameter) models 

We investigated whether estimating fewer parameters 
would lead to an increase in power. In some biological 
systems it may be reasonable to assume from the begin- 
ning that some processes are symmetrical. To examine 
this we compared the results of the previously 
described 5 vs. 6 parameter tests with the results of 
tests involving reduced models of 3 vs. 4 parameters. 
The reduced test compared a four-parameter model, 
where two of the rates were constrained to be equal in 
both character states, for example (X 0 , Xi, u 0 =Ui, 
<7oi=<7io)» and a three-parameter model where all three 
rates were constrained to be equal in both character 
states. These scenarios were done for two tip bias sce- 
narios, one with a small bias (3:1 character state ratio) 
and one with a greater degree of bias (7:1 character 



state ratio) as seen in Additional file 1: Table S4. We 
simulated sets of 500 trees with both levels of bias for 
each of the three parameter pairs as well as sets with- 
out bias for a six/five and four/three model compari- 
sons as seen in Additional file 1: Table S4. 

Estimating parameters in asymmetrical scenarios 

Using BiSSE to estimate rates of speciation, extinc- 
tion and character change may be illuminating not 
only to understand the degree of any asymmetries, 
but also to distinguish which potential factor (biased 
speciation, extinction, or character change) is respon- 
sible for an observed excess of species with a par- 
ticular state. Rate parameters for unconstrained (six 
parameter) models were tabulated under a best-case 
scenario representing a small degree of tip bias and a 
worst-case scenario that included a high degree of 
tip bias with tree sizes of 500 taxa. Parameters were 
estimated from the same 500 trees and respective 
characters that were used to calculate the BiSSE likeli- 
hood difference (Table 1). With these simulated cases 
we asked whether the parameter values were estimated 
well, with a special focus on whether a bias in one 
process (e.g., extinction) might be confounded with a 
bias in another process (e.g., character change). 

Implementation 

The computer software used in this study was a 
refined version of the package Diverse [25] described 
in Maddison et al. [1], with these refinements already 
implemented in Mesquite [26]. Two of these refine- 
ments are described here. The simulation module 
("Evolving Binary Speciation/Extinction Character") was 
modified to generate trees more efficiently by means 
of a continuous approximation. The updated module 
calculates the rate of events on the tree as the product 
of the number of terminal branches on the tree and 
sum of rates for each event type. The time to the next 
event was drawn from an exponential distribution and 
the type and location (branch) of the new event were 
drawn from appropriately weighted uniform distribu- 
tions. Following this, all terminal branches were 
extended to the time point of the generated event. 
This process continued until the tree reached a limit- 
ing number of tips, or the unlikely event that all ter- 
minal branches became extinct. The second refinement 
is an enhanced parameter estimator that uses a numer- 
ical integrator that implements the RKF45 method 
[28]. The RKF45 method improved on the RK4 
method, used in Maddison et al. [1] by adaptively 
adjusting the step-size used in the integration process. 
Our implementation specified a starting step-size and 
subsequent changes in step-size were limited a range 
of l/10x to lOx the original step-size. 
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Additional file 1: Table SI. Power of asymmetrical speciation rate 
simulations. Remaining parameters were symmetrical for each simulation 
(cjoi = 0.01, g 10 = 0.01, p 0 = 0.03, p, = 0.03). Power is plotted in Figure 1 A. 
The observed percent of terminal taxa with State 0 is indicated by the 
mean value from 500 simulations. Table S2. Power of simulations for 
character rate change. Remaining parameters were symmetrical for each 
simulation {p 0 = 0.03, p , = 0.03, X 0 = 0.1 , A , = 0.1). Power is plotted in 
Figure 1 B. The observed percent of terminal taxa with State 0 is indicated 
by the mean value from 500 simulations. Table S3. Power of 
asymmetrical extinction rate simulations. Remaining parameters were 
symmetrical for each simulation (q m = 0.01, q K = 0.01, A 0 = 0.1, A , = 0.1). 
Power is plotted in Figure 1C. The observed percent of terminal taxa with 
State 0 is indicated by the mean value from 500 simulations. Table S4. 
This table lists statistical power of the BiSSE model for 500 simulations 
containing 3:1 and 7:1 biases in terminal states for varying tree sizes for 
likelihood comparisons of power in four versus six parameter models. 
Power is plotted in Figure 2. Using the stationary frequency formula, in an 
iterative calculation, we obtained ratios of rates necessary to generate a 
low bias representative (3]'A 0 ) ar| d high bias representative (7y1o) tip ratios 
using values symmetrically placed around base rates (A = 0.1, u = 0.05, 
and q = 0.005). For the low bias, rate ratios were 1.1425, 1.3046 and 3.0 
for speciation, extinction and character change respectively, yielding 
simulated rates for speciation A o =0.0936, A,= 0.10689, for extinction, 
u 0 =0.04378, u ,= 0.05711, and for character change q 0 = 0.00289, 
c?,= 0.00866. For the high bias, rate ratios were 1.407, 1.960, and 7.000, 
yielding simulated rates for speciation Ao=0.0843, A^O.1186, for extinction 
u 0 =0.0357, u,= 0.07, and for character change q 0 =0.00189, q,= 0.01323. 
Simulated rates without bias were set to their base rates. 

Additional file 2: Figure SI. Parameter estimations of (a) speciation, 
(b) character change, and (c) extinction under different tree sizes and 
degrees of asymmetry in speciation rates with corresponding tip ratios. 
Point of intersection between red lines represents known values. 



